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ABSTRACT 



Investigation of two-dimensional transonic flows is 
extended to axi- symmetric problems. This is of considerable 
practical interest, for example, with regard to missiles or 
aircraft engines which approximate much more closely bodies of 
revolution than two-dimensional bodies. The main concern with 
axi- symmetric flows lies not only with the complexity of the 
governing nonlinear partial differential equation which is 
mixed of elliptic-hyperbolic type but also with the lack of a 
general method for accurately solving this type of equation. 
We solve the nonlinear transonic equation using separation 
of variables technique, which yields two nonlinear ordinary 
differential equations. The x-dependence can be integrated 
numerically, and the solution for the r-dependence can be 
obtained using the expansion method originated by Van Dyke. 
This works well with only three terms in the expansion. The 
sonic solution of these equations is obtained analytically 
since both equations are simple enough to be integrated for 
this case (1^=1. 0). The small parameter (1-M^ 2 ) plays an 
important role in specifying the shape of the boundary 
surfaces for external axi-symmetric steady flow of interest 
for design. A Navier- Stokes solver was used to compute the 
inviscid flow to confirm our results in the region over the 
surface where the small perturbations apply. 
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I. 



INTRODUCTION 



Transonic flow is a classical problem in gas dynamics and 
yet not many exact solutions are available. The main 
difficulty lies, of course, in the nonlinearity of the 
equations and in the boundary conditions . [Ref . 1] 

It is of interest to extend two-dimensional flow solutions 
to axi- symmetric problems, which are more practical for 
designing bodies of revolution such as missiles or aircraft 
engines. The goal is to design a shock- free body or surface so 
we can minimize the penalty of the wave drag that is caused by 
shock waves. The small perturbation technique may be 
implemented for simplifying the governing equations and the 
mathematical procedures for solving the resulting equations in 
nonlinear flows. This research shows a solution to nonlinear 
two-dimensional, external flow, over an afterbody axi- 
symmetric surface. Beginning with the transonic equation for 
an axi -symmetric body in Chapter II and applying separation of 
variables gives two nonlinear ordinary differential equations 
which lead into two exact solutions. One of the nonlinear 
ordinary differential equation can be integrated numerically, 
and other one can be solved using the outer expansion method 
by Van Dyke [Ref. 2], Chapter III shows the derivation of 
pressure coefficient and local Mach number for different 
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surfaces. Supersonic boundary surfaces are obtained 
implicitly in Chapter IV for M„= 1.05,1.1, and 1.2 and they are 

examined by a 3D Navier-Stokes solver, run as an Euler solver 
for this axi- symmetric problem in Computational Fluid Dynamics 
(CFD) . In addition, the sonic surface is derived analytically 
in Chapter V. This sonic boundary surface is also examined 
using 3D Navier-Stokes solver in CFD to verify the results. 
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II. 



TRANSONIC EQUATION FOR AXI - SYMMETRIC BODY 



The small perturbation, non-linear, axi- symmetric, 

transonic equation is written as [Ref. 3] 






1 _d_ 
i dr 



(r<p r ) 



XX 



( 1 ) 



Where the non-dimensional axial, x, and radial, r, velocity 
potentials have been defined as [Ref. 1] 

<P X = mL (y+l)-^, <p r = Ml (y+l) (2) 



where is free stream Mach number and y is ratio of heat 
capacities . 

An exact solution to Equation (1) , the transonic equation, 
can be found from 

<p (x,r) = <p s (x,r) + ( 1 -Ml ) x (3) 

Which patterned after the two dimensional case given by 
Biblarz [Ref. 1]. Here <p s (x,r) is a separable solution for 
the velocity potential. 
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The transonic equation in axi- symmetric form, Equation (1) , 
may be separated by letting the potential function <p (x,r) 
equal to 

<P (x,r) =Ux)( (r) + (1 -Ml)x (4) 



Substituting the above function into the transonic equation 
results in two ordinary, second order, non-linear differential 
equations 

(5) 

dx dx 2 

+ -XC 2 = 0 (6) 

dr 2 r dr 



where X is the separation constant. 

The solution to the first O.D.E. (5) is obtained by multiplying 
both sides by d\/dx, 



*L\ 

dx 



dj dH)_ dl 

dx dx 2 ) dx 



(X l) =0 



(7) 



or 



_d_ 

dx 



1 

3 





= 0 



(8) 
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Thus 



dx 




Rearranging 



dx = 



d5 



l ^ 2 



+ a 



1 




and integrating 




( 9 ) 



( 10 ) 



(ID 



where a and x a are integration constants. 

The solution to the second non-linear O.D.E. (6) , is obtained 
by an outer expansion method. [Ref. 2] 

C(r)=-i- + (1 -Ml)f x (r) + (1 -Ml) 2 f 2 {r) + . . . (12) 

Xr z 

Where (1-Mi) represents a small parameter and the first term 
is the purely sonic solution. 



t 
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By taking first and second derivatives and substituting 
them into the original 2nd order, non-linear, O.D.E.(6), the 
solution becomes 



C (r) = 



Ar : 



4 + 1 1 -Ml \aXr ( ^ +2) + HJhH . a 2 ( ^ +2) 



( 13 ) 



where a is constant. 



Equations 11 and 13 are solutions to the differential 
equations 5 and 6. Boundary conditions are needed to 
determine the constants a , c x , a, and X . 

The constants C lf a, and X can be shown to be related by the 
expression [Ref.l] 



' cA a 1 - 2 * 26 

k X ) X 0 ' 7574 



1 . 08X1CT 2 



( 14 ) 



and 



a = + Cjl-itfil 1 - 7574 



( 15 ) 



The positive sign above is for M„ a 1.0 and the negative sign 
is for < 1.0. 
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Upon inserting Eqn. (15) into Eqn. (11) , we have 




dt 



+ C, |l -Ml I 1 - 7574 



( 16 ) 



Factoring out C x , we obtain 



__i 

x-x 0 = Ci 3 






l-Mil 1 - 7574 ] 3 



( 17 ) 



Introducing new variables 



l ■ i 



3X 

12 C, J 
i/ 



( 18 ) 



(\/F + 2) = (\/^ + 2) ^ 



( 19 ) 



Equations (17) and (13) become 



x-x 0 




( 20 ) 
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c= (aX ) 0 - 4142 



Xr 



T" 2 



4 + [l-Mf ]f <^+ 2 > + (1 — r 2 '^ 21 
1 1 50.63 



( 21 ) 



Finally, defining two new variables 



3 X \ 2 






( 22 ) 



: - 



C >■ 

(a X ) 0 - 4142 



( 23 ) 



Equations 20 and 21 become 



x= 



dl 



[? 2 *|1-W. 2 | 1 - 7574 ]' 



where x 0 =0. 



( 24 ) 



r 2 



+ U-M 2 \ r ( v /F ^ 2) * — 2=21^2) 

1 50.63 



( 25 ) 



Equation (24) will be numerically integrated and plotted in 
Fig. (1) as l ( x ) vs. St , and Equation (25) will be evaluated 
and plotted in Fig. (2) as £ (r) vs. f for Af„= 1.05, 1.1, 

1.2. These results will be used later to obtain transonic 
surfaces with their corresponding C p and local Mach number 
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profiles. In Fig. (2) , Equation (25) is plotted up to its 
minimum value with increasing r. Thereafter, a constant 
value is patched in to the right of the + mark. This is done 
because of the need to keep £ constant reflecting the 
necessary behavior of the function for range r [Ref.l]. 
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III. 



PRESSURE COEFFICIENT AND LOCAL MACH NUMBER 



The pressure coefficient for a 
perturbations is given by [Ref. 4] 



flow with small 



2u y 
“ + 



(1 -Ml) 





(26) 



The linearized pressure approximation for axi- symmetric flow 
is 



a 



~ 2 U x 

U. 



(27) 



Recall the axial velocity potential, Equation (2), Chapter II, 



cp x = Ml (Y+l) 



(28) 



thus 



_^x 



<Px 

Ml (y + 1) 



(29) 
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substitute Equation (29) into Equation (27) , yields 



-2 <Px 
Mi (Y + l) 



( 30 ) 



We need to take the derivative of the potential function with 
respect to x in equation (4) . 



<P* 



- r 

^ dx 



+ (i - Mi) 



( 31 ) 



Rewriting equation (9) with the constant a inserted becomes 



d£ 

dx 



3_X 

2 



V + q | (i-*d) I 






1.7574 



( 32 ) 



Recall Equation (18) 



l = 5 




( 33 ) 



Inserting Eqn. (33) into Eqn. (32) and factoring C 1 out, yields 



dx 



.1 j. 

Q 3 [ £ 2 + | (1 -Mi) I 1 - 7574 ] 3 



( 34 ) 
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Recall Equation (23) and rearrange 



c = (a X ) 0 - 4142 ? 



( 35 ) 



Rewriting Equation (31) the potential function as 



<Px = 



_ (ak) °- 4142 , — 3 



£ C x 3 [ £ 2 + | (l-Mf ) I 1 - 7574 ] 3 + {1-mI) (36) 



Thus 



<P* = 



, 0 . 4142^3 



0.5858 



l [ Md -id) l 1 - 7574 ] 3 ^!-^ 2 ) 



( 37 ) 



Recall Equation (14) Chapter II 



£i 

k 



a 1 - 2426 
X 0,7574 



1 . 08xl0' 2 



( 38 ) 



then finally Equation (37) becomes 



<p x = 0.2208 £ [ l 2 + | (1-Af.) I 1 - 7574 ] 3 + (1 -m2) 



( 39 ) 



Rewriting the coefficient of pressure 



C = 



-2 



P Ml{ y+1) 



0.2208 £[£ 2 + | (1 -Ml) I 1 - 7574 ] 3 + (1 -Ml) 



( 40 ) 
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The previous expression for C p will be used further in 

determining the local Mach number. 

The perturbation velocities are assumed to be vary small 
compared to the undisturbed uniform velocity U„. [Ref. 5] 

Hence at a point near the body, the velocity vector V for 
2-D flow is given by 

v = i {U m + u x ) + j u r (41) 



We can obtain the local Mach number M in terms of the 
perturbation velocities and the local speed of sound c as 

M 2 = Jiii = + U *> 2 + U r (42) 

C 2 C 2 



Neglecting the higher order ratios of the free stream to 
perturbation velocities, since they are considerably smaller 
than unity, we have 



M 2 



u: 



1 + - 



2U> 



(43) 
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The ratio of the local speed of sound to the undisturbed 
uniform speed of sound becomes 



c 



2 




( 44 ) 



or 



c 2 = TjT m 
cl ’ T 0 /T 



( 45 ) 



Where T a is the stagnation temperature. 

We can write temperature ratios as functions of Mach numbers 
as 




1 + ( y - 1) Ml 
2 

1 + (Y-D M 2 

2 



( 46 ) 



Substituting Equation (45) into Equation (46) and rearrange, 
we have 



M 2 = Ml 



1 + 
^ 



£^Vi + Jyz1> 

u. Jl 2 




Y-l 

2 





( 47 ) 
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from Equation (27) 



2 

U m 



Rearranging and solving for local Mach number 



M 2 



MZ (1 -C p ) 



1 + — — 

2 



Ml C„ 



This equation will be used to show the local Mach 
each surface at M„ = 1.05, 1.1 and 1.2 [Ref. 6]. 



( 48 ) 



(49) 



number for 
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IV. 



SUPERSONIC BOUNDARY SURFACES 



The boundary conditions which need to be applied require 
that the gradient of <p vanish far ahead of the body and that 
the flow be tangential to the surface [ Ref. 7] . In terms of 
perturbation velocities the boundary conditions become 




( 50 ) 



Recall the modified velocity perturbation potential Eqn. (2) 



<P r = Ml ( Y + 1 ) 



( 51 ) 



then, we have 



= <Pr 

u - Ml (y + 1) 



( 52 ) 



by substituting Equation (52) into Equation (50), we obtain 



I dr\ = <Pr 
\ dxj Ml (y + 1) 



( 53 ) 
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By differentiating Equation (4) with respect to r gives 



<Pr 




( 54 ) 



Recalling Eqn. (18) and substituting Eqn. (21) into Eqn.(23), 
and taking the derivative with respect to r, the above 
equation becomes 



<Pr 



0 . 0848 




( 55 ) 



Substituting Equation (55) and (13) into equation (53), we 
have 




°- 0848 l + 2.8284 |l-A*f | f 1 ' 8284 

Ml{ y + D 

+ 0 . 1512 (1 -Ml) 2 i 6,657 



( 56 ) 



It can be noticed here that l and £ are given implicitly as 
function of x and f respectively. 

From Eqns. (19) and (22), we can rearrange Equation (56) as 



dr 

dx 



3 .26x10 2 £ + 2.8284 f 1,8284 

M 2 (y+i) Lr 3 



+ 0.1512 (1-Mi) 2 f 6,657 



( 57 ) 
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Further arrangment of the previous equation yields 



dr 

-2 . 82 84 1 1 | f 1 - 8284 - 0 . 1512 (1-^f ) 2 f 6 - 657 

- 3 - ' 26 ^ 0 ' 2 l dx 
mI{ y+i) 



( 58 ) 



Now, we can compute the exact boundary surfaces out of the 
equation above. Starting with the left hand side of the 
equation 



Z(f) = 



8-2 . 8284 \l-Ml\r 



f 3 

_ M 2| ~ 4 .8284_ 



0 . 1512 (1 -MZ) 2 f 



-M 2 \ 24=9.657 



( 59 ) 



By defining equation (59) we can plot Z (r) vs. r in Fig. 

(3) for = 1.05, 1.1, and 1.2. 

The minimum values of r a for different Mach numbers can be 
expressed as r a = r Q (AfJ which has been found to be [Ref. 1] 



1.2074 

| 1-^|°- 2071 



( 60 ) 
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Fig. (4) represents r Q vs. M „ and shows a symmetry close to 
Af«,=1.0 where r Q =°°. So for our purpose to represent the 
condition sufficiently close to = 1.0, we will choose a 

finite value of r 0 (perhaps as 4.0). 

Defining the right hand side of Equation (58) as 

Kix) = zl^xicrl f j ^ (61) 

«-< Y+l) J„ 



taking the derivative of Equation (24) and substitute it into 
Equation (61) , yields 



K(x) = 



- 3 . 26x1 0 ~ 2 
Ml (y + l) 



f. 



Idl 



• [P-Hi-M-I 1 ' 7574 



] 3 



(62) 



integrating by parts, we obtain 



K(x) = 



-2 . 45x1 O' 2 



(y + l) Ml 



(V - | !■ -Ml I 1 - 7574 ) 3 + 1 1 - -Ml I 1 - 1716 



(63) 



the equation above will allow us to plot K(x)vs.x in 
Fig. (5) for Af„ = 1.05, 1.1, and 1.2. Where K (x) < 0 

for 1.0 and K (x) > 0 for M„ <1.0. 
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Therefore, integrating both sides of equation (58) becomes 




I 




X 




mI{ y+i) 



Equation (64) is numerically integrated and which then allows 
us to determine the boundary surfaces in dimensionalize and 
non-dimensionalize (normalize) form as shown in Figs. 6 and 7 
for M„ = 1.05, 1.1, and 1.2. Now, we can determine C p from 

Eqn. (40) and local Mach number from Eqn. (49) for the three 
surfaces obtained at different Mach numbers. 

In Figs. 8-13, the supersonic boundary surfaces for 
M 00 =l . 05 , 1 . 10 , and 1.20 are plotted with their pressure 
coefficients and local Mach profiles in dimensionalized and 
non-dimensionalized (normalized) form. These results will be 
verified later using a 3-D Navier-Stokes solver in 
Computational Fluid Dynamics (CFD) . 
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V. 



SONIC BOUNDARY SURFACE 



The sonic flow (Af„=1.0) can be derived in explicit form 

from Equations (24) and (25) . The resulting equations for 
sonic flow are 



x = 




( 65 ) 



and 




( 66 ) 



Rearranging Equation (65) 



x = 




2 

3 dz 



( 67 ) 



since Equation (67) is integrable, we obtain 

x = 3 V 
or 




( 68 ) 



( 69 ) 
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We can write Equations (11) and (13) in sonic form as 



x = 




dt, 




i 

3 



and 



C = 



4 

X r 2 



integrating Equation (70) , we have 

Z(x) 



( 70 ) 



(71) 



(72) 



From Equation (4) , we can write the sonic perturbation 
potential as 

<p (x, r) = Z (x) C (r) (73) 

Substituting Equations (71) and (72) into Equation (73) yields 

<pU,r) = (74) 

9 r 2 

Equation (57) can be written for Af„=l . 0 as 

dr _ - 0 . 2608 \ lncs 

~dx = — fit? ( 1 
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Rearranging Equation (75) , we have 



r 3 dr =-0 . 1087 l dx 



( 76 ) 



Therefore, integrating both sides of Equation (76) , we obtain 



f x 3 dx = - 0 . 1087 f l dx 

Jo 



( 77 ) 



Inserting Equation (69) into Equation (77) yields 



f 4 -f 0 4 = - 0.0040 x 4 



( 78 ) 



where f 0 = 4.0 has been chosen as an asymptotic value to 

represent the condition sufficiently well at Af„=1.0. Equation 

(78) will allow us to determine sonic boundary surface. 

The pressure coefficient can be written from Eqn. (40) as 



C„ = 



_ ~ 0 . 4416 ^ 3 



(Y + l) 



( 79 ) 



in case of M„=l . 0 . 
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The local Mach number can be written from Equation (49) for 
the sonic surface as 



M 2 = 



t 1 - C ,) 



+ 



v-1 

2 




( 80 ) 



The sonic surface is represented dimensionally in Fig. (14) 
with the pressure coefficient and local Mach number profiles, 
and normalized in Fig. (15) with r 0 =4 . 0 for (M„=1.0). The 

sonic boundary surfaces, in both forms, are plotted with the 
supersonic boundary surfaces in Figs. 16 and 17. 
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VI. 



COMPUTATIONAL FLUID DYNAMICS (CFD) 



The history of Computational Fluid Dynamics (CFD) is 
closely tied to the rapid advances in digital computer which 
has a great impact on problems of design in modern engineering 
practice. These problems can now be solved at very little cost 
in a few seconds of computer time which would have taken years 
to work out with the computational methods and computers 
available twenty years ago. [Ref. 8] 

The CFD will be used to compute the axi- symmetric flow 
over the boundary surfaces obtained by the small perturbation 
method. 

The computer program GRAPE [Ref. 9], an acronym derived 
from Grids about Airfoils using Poisson's Equation, has been 
written by R. Sorenson at Ames Research Center to generate 
two-dimensional finite difference grids about airfoils and 
other shapes by the use of the Poisson differential equation. 
Outer and inner boundaries are specified for the C-type grid, 
where the surface of the body is treated as the inner 
boundary. Important characteristics in grid generation are 
the ability to specify the spacing between mesh points at the 
boundary, in the direction normal to the boundary, and the 
control of the angles with which mesh lines intersect the 
boundaries which is known as orthogonality. 
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The grid size for the sonic surface (M^l.O) is 107 x 60 and 
for the supersonic surfaces (M^-l.l, 1.2) is 115 x 60 as 
shown in Fig. (18) . Since the surfaces are symmetrical, half of 
the grid or eventually the lower surface was rotated for 11 
planes plus two more for symmetry to generate an axi- 
symmetric after body surface. 

The OVERFLOW program [Ref. 10] developed by Ames Research 
Center which uses 3-D Navier- Stokes and Euler solver for 
viscous/ inviscid flow was applied. The results are shown in 
Figs. 19-21 for three different boundary surfaces at = 1.0, 
1.1, and 1.2 with 1500 iterations. It was found that the 
density residuals decreased by more than one order of 
magnitude over 1500 iterations. After 3000 iterations, the 
solution converged by two orders of magnitude. 

In case of the sonic surface Fig. (19) with M^l.0, the 
Mach lines contour shows the shock is forming downstream at a 
local Mach M=1.65 at which the flow becomes subsonic 
downstream of the shock. This result complies with the small 
perturbation transonic solution obtained earlier in Fig. (15) . 
In other words, the flow is shock free over the sonic surface 
in for M^l . 0 . 

For the supersonic boundary surfaces M^l. 10, and 1.20, we 
can see that the flow starts at M=1.10 and M=1.20 respectively 
and follows the afterbody surface until it shocks. The shock 
is formed at local Mach M=1.75 and 1.85 respectively as shown 
in Figs. 20 and 21. 
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On the other hand, these results agree with the small 
perturbation solution over the transonic range as plotted in 
Figs. 11 and 13 where the local Mach is represented by the 
dash line and it is increasing towards the afterbody surface 
until it shocks. This indicates that the flow is shockless 
over these supersonic boundary surfaces in the transonic 
regime. 

The velocity vectors are plotted in Fig. (22) for M co =1.0 
which shows that the flow is following the afterbody surface 
until the shock, and then flow separation takes place due to 
the steep pressure rise and the introduction of numerical 
viscosity into the flowfield at this location. 

These boundary surfaces obtained are of interest for 
design of body of revolution such as missiles afterbody or 
aircraft engines, in which a patching technique may be applied 
to these surfaces to get the best shockless surface in the 
transonic range. 
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VII. 



CONCLUSION AND RECOMMENDATION 



We have shown in this research a solution to the non- 
linear transonic small perturbation equation by implementing 
the separation of variables technique. The x- dependence was 
integrated numerically and the r- dependence was solved by an 
outer expansion method. It was found that the parameter (1- 
M^ 2 ) has strong effect in specifying the shape of the boundary 
surfaces of interest for design. A 3-D viscous/ inviscid 
Navier- Stokes/ Euler solver in Computational Fluid Dynamics 
(CFD) confirmed our results obtained from the small 
perturbation technique for the boundary surfaces at M^l.O, 
1.10, and 1.2 0 . In other words, we achieved our goal of 
designing shockless surface in the transonic range. 

This research may be extended using a patching technique 
to obtain the best shockless surface in the transonic flows. 
Patching criteria will have to be developed based on the 
mathematics and the flow constraints. Also, connection between 
small perturbation solution and CFD might eventually be done 
internally in the computer. The small perturbation is the most 
efficient method to define or to design axi- symmetric 
afterbody transonic surfaces, however CFD can be used to 
predict the performance of these aerodynamic surfaces. 
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APPENDIX A 
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Figure 1. Numerical integration of Eqn . 24. 
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Figure 2. Numerical solution of Eqn. 25. 
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Figure 3. Numerical solution of Eqn. 59. 




33 



Figure 4 . Numerical solution of Eqn. 60. 
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Figure 5. Numerical solution of Eqn. 63. 
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Figure 6. Numerical solution of Eqn. 64 
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Figure 7. Numerical solution of Eqn. 64 (normalized). 
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Figure 8. Cp and local Mach for 



Cp & local Mach FOR M <x> = 1.05 , ro = 1.935 
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Figure 9. Cp and local Mach for H. = 1.05 (normalized) 
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Figure 10. Cp and local Mach for M. = l.io. 
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Figure ll. Cp and local Mach for M, = l.io (normalized) 
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Figure 13. Cp and local Mach for M,, = 1.20 (normalized). 
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Figure 14. Cp and local Mach for 
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Figure 15. Cp and local Mach for M,, = 1.0 (normalized) 
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Figure 16. Boundary surfaces for 
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Figure 17 . Boundary surfaces for M. = 1.0/1.05/1.1,1.2 (normalized) . 
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Figure 18. Afterbody axi- symmetric grid. 
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FIGURE 2D. HACH CONTOURS FOR M 
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FIGURE 21. MACH CONTOURS FOR 
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Figure 22. Velocity vectors for 



APPENDIX B 
COMPUTER PROGRAMS 



★ * 

* This program is designed to calculate KSI(X), * 

* K(X), and X using numerical integration based * 

* on trapezoidal rule to solve Eqn. (24) and * 

* plotting the output as shown in Figs.l & 3 . * 

* * 
**************************************************************** 

PROGRAM KSI 

REAL M(3) ,X( 0:401, 3) , A, Al , B (401 , 3 ) 

+ , KS I , P , FUNC , H , DD ,Q,L,K(401,3) ,N,N1 

INTEGER YY 

OPEN (UNIT=9 , FILE= ' KSI ' , STATUS= ' UNKNOWN' ) 

PRINT *, 'ENTER LOWER & UPPER BOUND, # OF INTERVALS, 

+#OF DATA SET' 

READ *, A, Al, N, N1 

DO 10 1=1,3 

PRINT * , ' ENTER MACH NO . ' 

READ *, M ( I) 

M ( I) = 0. 9 + 1*0. 1 

IF (M(I) .LT.1.0) THEN 
P=-l. 

ELSE 
P=l. 

ENDIF 

DO 20 J=1 , N1 
B ( J, I) = J*A1/N1 
H = (B ( J , I) -A) /N 
AREA = 0. 

K( J, I) = (-0.0102/ (M(I) **2) ) * ( (ABS ( (B ( J, I) **2) - 
+( (ABS (l-M(I) **2 ) ) **1.7574) ) ) ** (2./3 . ) + 

+ (ABS (l-M(I) **2) ) **1.1716) 

DO 30 L=1 , N 
KSI = A + (L-0 .5) *H 

DD = (KSI**2+P* (ABS (l-M(I) **2) ) **1.7574) 

IF (DD.GT.0.) THEN 
FUNC = 1./ (DD** (1./3 . ) ) 

Q = 1. 

ELSE 



100 

C 
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30 

20 

10 

50 

40 

200 

110 



FUNC = 1./ ( (-1*DD) ** (1./3. ) ) 

Q = -1. 

ENDIF 

AREA = AREA + H*FUNC*Q 
X(J,I) = AREA 

CONTINUE 

CONTINUE 

CONTINUE 

DO 40 J=1 , N1 

PRINT 50, (B(J,I) ,X(J,I) ,K(J,I) ,1=1,3) 

FORMAT (IX, 3 ( 1F6 . 3 , 2F10 . 5 , IX) ) 

CONTINUE 

PRINT* , ' TRY AGAIN ? ENTER 1 FOR YES , 2 FOR DATA 
+FILE, OTHERS FOR NO' 

READ *, YY 
IF (YY.EQ.l) THEN 
GO TO 100 
ELSE 

IF (YY.EQ.2) THEN 
DO 110 J=1 , N1 

WRITE (9,50) (B(J,I) ,X(J,I) ,K(J,I) ,1 = 1,3) 

CONTINUE 
GO TO 200 
ENDIF 
ENDIF 
END 
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★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★a-* 



★ ★ 

* This program is written to calculate ZETA(r) and * 

* r using Eqn. (25) and plotting the output in Fig. 3 . * 

* * 



★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★•A: 



program zeta 

real M(3 ) , zeta (0 : 14 , 3 ) , r (14 ) 

open (unit=9 , f ile=' zz ' , status=' unknown' ) 

do 10 1=1,3 

c M(I) =0. 7+1*0. 1 

print *, 'enter mach no.' 
read *, M(I) 

do 20 J=1 , 14 
zeta (0,1) =1000 
print *, 'enter r values' 
read *, r(J) 
c r ( J) =0 . 1*J 

zeta (J, I) = (4/ (r (J) **2) ) + (abs (l- (M(I) **2) ) * 

+ (r (J) **2. 8284) ) + ( (1- (M(I) **2) ) **2) * (r (J) **7 . 657) /50 . 63 
if (zeta (J, I) .GE. zeta (J-l, I) ) zeta ( J, I) =zeta ( J- l , I) 

20 continue 

10 continue 

print 30, M 

30 f ormat (llx, 5F10 . 4 ) 

do 40 J=1 , 14 

c print 50, r ( J) , ( zeta ( J, I ) , 1=1 , 4 ) 

write (9,50) (r ( J) , zeta ( J, I) , 1=1 , 3 ) 

50 format (4x, 3 (2F10 . 4 , 2x) ) 

40 continue 

end 
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★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it*** 



★ * 

* This program is written to calculate Z(r) and r * 

* using Eqn. (59) and showing the output in Fig. 4 . * 

* * 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★it** 



PROGRAM ZR 

REAL M(4) , R ( 1000 ) ,Z(0:1000,4) 

OPEN (UNIT=99 , FILE= ' ZR1 ' , STATUS =' UNKNOWN' ) 

DO 10 1=1,4 

C M(I) = 0.7 + 1*0.1 

PRINT *, 'ENTER MACH NO. ' 

READ *, M ( I ) 

DO 20 J=1,1000 
Z(0,I) = 0.0 
R ( J) = 0.001*J 

Z(J,I) = (R(J) **3) / (8 

(2.8284* ( R ( J) **4.8284) *ABS (1- (M(I) **2) ) ) - 

+ (0 . 1512* (R (J) **9 . 657) * ( (l-M(I) **2) **2) ) ) 

IF (Z(J,I) .LE.0.0) Z (J, I) =Z (J-l, I) 

20 CONTINUE 

10 CONTINUE 



DO 40 J=l,1000 

C PRINT 50, R ( J) , (Z(J,I) ,1=1,4) 

WRITE (99,50) R(J), (Z(J,I), 1=1,4) 

50 FORMAT ( 4X, 5F10 . 6 ) 

40 CONTINUE 

END 
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This program is written to determine Rmin for 
transonic Mach numbers using Eqn. (60) and the 
output is shown in Fig. 5 . 



**************************************************************** 



PROGRAM RMIN 
REAL M,R 

OPEN (UNIT=88 , FILE-' RMIN' , STATUS =' unknown' ) 

DO 100 M = 0.8,1.2,0.0001 
IF (M. LT .1.0) THEN 

R=1 .2074/ ( (1-M**2) **0.2071) 

ELSE 

IF (M. GT .1.0) THEN 

R=1 .2074/ ( (M**2-l) **0.2071) 

ELSE 

R = 7.0 



★ 



★ 



★ 



C 



80 

100 



ENDIF 

ENDIF 

PRINT 80, M, R 
WRITE (88,80) M, R 
FORMAT (IX, F7 . 5 , 7X, 1F10.6) 
CONTINUE 
END 
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**************************************************************** 



* * 

* This program is designed to numerically integrate * 

* Z(r) ,Eqn.(64), based on trapezoidal rule . * 

* * 
**************************************************************** 



PROGRAM IZR 

REAL M(3) , R (401 , 3 ) , A, Al, IZ (0 : 401, 3) , H, FUNC, N, N1 , Ro 
INTEGER YY 

OPEN (UNIT=55, FILE-' IZR' , STATUS =' UNKNOWN' ) 

100 DO 10 1=1,3 
C M ( I ) = 0.7 + 1*0.1 

PRINT *, 'ENTER MACH NO. ' 

READ *, M ( I ) 

PRINT *, 'DETERMINE THE VALUES OF MACH #',M(I) 

PRINT *, 'ENTER Rmin ,A' 

READ *, A 

PRINT *, 'ENTER THE UPPER LIMIT OF R , Al ' 

READ *, Al 

PRINT *, 'ENTER # OF INTERVALS , N' 

READ * , N 

PRINT *, 'HOW MANY DATA SET DERIVED ? Nl' 

READ *, Nl 

PRINT *, 'INPUT Ro FOR THIS MACH # , Ro' 

READ *, Ro 

DO 20 J=1 , Nl 
R ( J, I) = J*A1/N1 
H = (R ( J, I) -A) /N 
AREA = 0. 

DO 30 L=1 , N 
RR = A+ (L-0.5) *H 

FUNC= (RR**3 ) / (8- (2 . 8284* (RR**4. 8284) *ABS (1- (M(I) **2) ) ) 
+ - (0.1512* (RR**9 .657) * ( (1-M (I) **2 . ) **2 . ) ) ) 

AREA = AREA + H*FUNC 
IZ(J,I) = Ro - AREA 

30 CONTINUE 

20 CONTINUE 

10 CONTINUE 

DO 40 J=1 , Nl 

C PRINT 50, (R(J,I) ,IZ(J,I) ,1=1,3) 
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WRITE (55, 50) (R ( J, I) , IZ ( J, I) , 1=1,3) 

50 FORMAT (IX, 3 (F6 . 3 , F9 . 6 , 3X) ) 

40 CONTINUE 

PRINT * , ' TRY AGAIN ? ENTER 1 FOR #YES#, OTHERS FOR #NO# ' 
READ *, YY 

IF (YY.EQ.l) GO TO 100 
END 
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**★*****★******★***★*******★★★★★★*★★★★*★*★★**★**★★***★***★***★** 
* ★ 

* This program is using numerical integration to 

* determine the supersonic boundary surfaces , 

* Egn. (64) , and calculating pressure coefficient 

* and local Mach number profiles, using Eqns. (40 

* & 49) . The output is shown in Figs. (6-13) . 

* 

**************************************************************** 
PROGRAM CPP 

REAL X(401,9) , Y ( 401 , 6 ) , XR (401) , YR (401) ,M,MX, XYMAT 
INTEGER I , J, K, L, IWRITE 

0PEN(UNIT=52,FILE=' cpl5.dat' , STATUS =' UNKNOWN' ) 

OPEN (UNIT-53, FILE-' cpll.dat' , STATUS =' UNKNOWN' ) 
OPEN(UNIT=54,FILE=' cpl2.dat' , STATUS =' UNKNOWN' ) 

OPEN (UNIT=77, FILE-' cpp.dat' , STATUS- ' UNKNOWN' ) 

OPEN (UNIT-78, FILE-' ml05 .dat' , STATUS- ' UNKNOWN' ) 

OPEN (UNIT-79, FILE- 'tirel.dat' , STATUS- ' UNKNOWN' ) 

OPEN (UNIT-80, FILE-' tire2 .dat' , STATUS- ' UNKNOWN' ) 

100 OPEN (UNIT-1, FILE- 'KSI. DAT' , STATUS- ' OLD ' ) 

OPEN (UNIT-2, FILE-' IZR. DAT' , STATUS- ' OLD' ) 

REWIND ( UNIT- 1) 

REWIND (UNIT-2) 

IWRITE =78 
MACP-52 

READ (1,*, END-50) ( (X ( I , J) , J=1 , 9 ) , 1=1 , 401 ) 

50 READ (2,*, END-60) ( (Y(K,L) ,L=1,6) ,K= 1,401) 

60 DO 10 J=2 ,6,2 

PRINT *,' Ro is rmin ; EPSILON is the accuracy ' 

PRINT *, 'ENTER MACH NO., Ro , EPSILON' 

READ *, M, RO , EPS 
WRITE (77,*) M, RO , EPS 
XM = 1 -M**2 
L = J/2*3 

DO 30 I = 1,401 
DO 20 K = 1,401 

IF(Y(I, J) .LT. 0.001) THEN 
GO TO 10 
ELSE 

XYMAT = ABS ( Y ( I , J) - ABS (X (K, L) ) ) 

IF (XYMAT.GT.EPS) THEN 
GO TO 20 
ELSE 



* 

* 

* 

* 

★ 

* 
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XR(K) = X (K, L- 1) /RO 
YR ( I ) = Y ( I , J- 1 ) /RO 

ZETA = (4 . / (Y ( I , J- 1 ) **2) ) + ABS (XM) *Y(I, J-l) **2 
+ + XM**2* (Y(I, J-l) **7.657) /50. 63 

Cp = ( - 2 / (M**2 *2 . 4 ) ) * ( 0 . 2208*ZETA* (X (K , L - 2 ) **2 
+ (ABS (XM) ) **1.7574) **(l./3.)+XM) 

MX = SQRT (ABS ( (M**2* (1-Cp) / (1+0 . 2*M**2*Cp) ) ) ) 

C PRINT 90, X(K,L-1) ,Y(I,J-1) , Cp , MX 

WRITE (77, 40) X(K,L-1),Y(I,J-1) ,Cp,MX 
WRITE (MACP, 40) XR(K) ,YR(I) , Cp ,MX 
WRITE (IWRITE, 70) XR(K) , YR(I) 

GO TO 30 
ENDIF 
ENDIF 

20 CONTINUE 

30 CONTINUE 

IWRITE= IWRITE +1 
MACP =MACP+1 
10 CONTINUE 

40 FORMAT ( IX, 4 (F15 . 7 , 2X) ) 

70 FORMAT (1X,2F17.9) 

90 FORMAT ( IX, 4 (F15 . 7, 2X) ) 

PRINT * , ' TRY AGAIN ? 1 FOR YES ! , 2 FOR NO ! ! ' 

READ * , IY 

IF(IY.EQ.l) GO TO 100 
END 



.8284 

+ 
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**************************************************************** 
* * 

* This program is written to determine the sonic * 

* boundary surface M = 1.0, ro = 4.0 using Eqn. (78) , * 

* and calculating pressure coefficient and local * 

* Mach profile using Eqns.(79 & 80). The output is * 

* shown in Figs. 16 & 17 . Some output files are used * 

* in CFD . * 

* * 
**************************************************************** 



program sonic 

real M(0:160) , r ( 0 : 160) ,x(0:160) , zeta (0:160) ,ksi(0:160) , 
+ z (0:160) ,k(0:160) ,Cp(0:160) , RR ( 0 : 16 0 ) , YY ( 0 : 160 ) , 

+ xx ( 0 : 160) , xxx (0 : 160) 

open (unit=ll , f ile= ' sonicl . dat ' , status= ' unknown' ) 
open (unit=22 , f ile= ' sonic2 . dat' , status= ' unknown' ) 
open (unit=33 , f ile= ' sonic3 . dat ' , status= ' unknown' ) 
open (unit=66 , file=' sonic4 .dat' , status= ' unknown' ) 



c 


sonicl . dat 


= > x , r 


, Cp , M 


c 


sonic2 . dat 


= > x , r 


(horintal) 


c 


sonic3 . dat 


= > x , r 


(vertical) 


c 


sonic4 . dat 


=> x/ro , 


r/ro , Cp , M 



do 10 I - 0, 159 
x(I) = 1*0.1 
ro = 4.0 

r ( I) = (abs (ro**4 . - 0 . 004* (x (I) **4 . ) ) ) **0 . 25 

RR ( I ) = r(I)/ro 

YY ( I ) = - 1.0* RR ( I ) 

xx ( I ) = x(I) /ro 

xxx ( I ) = 3.975 - xx ( X ) 

zeta(I) = 4 . / (r (I) **2 . ) 

ksi(I) = (x(I) **3. ) /27.0 

z(I) = r (I) **3 . 

k ( I) = - 0.00101* (x(I) **4 . ) 

Cp(I) = - 0.1840 * zeta(I) * (ksi (I) **0.6667) 

M ( I ) = (abs ( (l-Cp(I) )/ (1 + 0.2*Cp(I) ) ) ) **0.50 
c print *, xx(I) ,RR(I) ,Cp(I) , M(I) 

10 continue 

do 20 I =0, 159 

write(ll,50) xx(I) , rr(I) , Cp(I) , M(I) 
write (22,60) xxx(I) , RR(I) , YY(I) 
write (33,70) xxx(I),YY(I) 
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50 

60 

70 

20 

80 

30 



write(66,5C) x(I) ,r(I) ,Cp(I) ,M(I) 
format (2x, If 8 . 4 , 2x, 3f 14 . 7 , 2x) 
format (2x, If 8 . 3 , 2x, 2f 12 . 6 , 2x) 
format (lOx, If 12 . 6 , 5x, f 12 . 6 ) 
continue 

write (34,80) (xxx (I) ,i= 0,159,3) 
write (34,80) (xxx(i) , i=156, 0, -3) 
write (34,80) (yy(i) ,i=0,159,3) 
write (34,80) (rr (i) , i=156, 0, -3) 
format (5 (lx, f 12 . 6, ' , ' ) ) 

do 30 I =158, 0 , -1 
write (33,70) xxx(I),RR(I) 
continue 
end 
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